\[\begin{align} S' & = \mu N^\star - \lambda S - \mu S \\ E' & = \lambda S - \kappa E - \mu E \\ I_S' & = p \kappa E - (\gamma_S + \mu_{I_S} + \delta_H) I_S - \mu I_S \\ I_A' &= (1 - p) \kappa E - \gamma_A I_A - \mu I_A \\ H' &= \delta_H I_S - (\gamma_H + \mu_H) H - \mu H \\ R' & = \gamma_S I_S + \gamma_A I_A + \gamma_H H - \mu R \\ D' &= \mu_{I_S} I_S + \mu_H H \\ \lambda &:= \frac{\beta_A I_A + \beta_S I_S}{N^{\star}} \\ N^{\star}(t) &= S + E + I_S + I_A + H + R \end{align}\]
| state | Stan | Initial Condition | Stan | Postulated |
|---|---|---|---|---|
| \(S\) | y[1] | \(S_0\) | y_init[1] | \(N - (I_{S_0} + I_{A_0} + E_0)\) |
| \(E\) | y[2] | \(E_0\) | y_init[2] | Unif(74, 2100) |
| \(I_S\) | y[3] | \(I_{s_0}\) | y_init[3] | 22 (cdmx data) |
| \(I_A\) | y[4] | \(I_{a_0}\) | y_init[4] | Unif(74, 2100) |
| \(H\) | y[5] | \(H_0\) | y_init[5] | 0 |
| \(R\) | y[6] | \(R_0\) | y_init[6] | 0 |
| \(D\) | y[7] | \(D_0\) | y_init[7] | 0 |
| \(I_s^{[cumulative]}\) | y[8] | \(I_{s_0}^{[cumulative]}\) | y_init[8] | 74 (cdmx data) |
\[\begin{align} Y_t & \sim \text{Poisson}(\lambda_t), \qquad \\ \lambda_t & = \int_0^t \rho \delta_E E \\ & p \sim \text{Uniform(0.3 0.8)} \\ & \kappa: \sim \text{Gamma(10,50)} \end{align}\]
| stan notation | symbol | Prior | Fixed |
|---|---|---|---|
| theta[1] | \(\beta_S\) | Normal(1, 0.3) | |
| theta[2] | \(\beta_A\) | Normal(1, 0.3) | |
| theta[3] | \(\kappa\) | Gamma(10, 100) | |
| theta[4] | \(p\) | Unif(0.3, 0.8) | |
| theta[5] | \(\delta_H\) | Gamma(10, 40) | |
| theta[6] | \(\mu_{I_S}\) | Gamma(10, 80) | |
| theta[7] | \(\mu_H\) | Gamma(10, 350) | |
| theta[8] | \(\gamma_S\) | Gamma(10, 350) | |
| theta[9] | \(\gamma_A\) | Gamma(10, 350) | |
| theta[10] | \(\gamma_H\) | Gamma(10, 350) | |
| y_init[2] | \(E_0\) | Unif(74, 2100) | |
| y_init[3] | \(I_S\) | ———– | 22 |
| y_init[4] | \(I_A\) | Unif(74, 2100) |
knitr::opts_chunk$set(echo = FALSE)
library(ggplot2)
library(plotly)##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(rstan)## Loading required package: StanHeaders
## rstan (Version 2.21.5, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
mu_ <- function(lambda, ylabel= "mu_h"){
# We suppose mu_. as a exponential r.v.
x <- seq(from = 0, to = 1, by=.001)
mu_h_ <- dexp(x, rate = lambda)
df <- data.frame(x=x, mu_h = mu_h_)
mean_ <- 1/lambda
mean_time <- 1/mean_
text <- paste(mean_time,"days")
p <- ggplot(data=df,
aes(x = x,
y = mu_h)) +
xlab("x")+
ylab(ylabel) +
geom_line(colour = "#000000") +
geom_area(alpha = 0.6, fill = "lightgray") +
geom_vline(xintercept = mean_, linetype="dotted",
color = "red", size=0.5)+
annotate(geom = "text", x = 1.30 * mean_, y=0, label = text,
color="black")
return(p)
}
#
gamma_ <- function(alpha, beta, ylabel="gamma_S"){
# We suppose mu_. as a exponential r.v.
theta <- 1 / beta
x <- seq(from = 0, to = 1, by=.001)
gamma_s_ <- dgamma(x, shape = alpha,
scale= theta)
df <- data.frame(x=x, gamma_s = gamma_s_)
mean_ <- alpha * theta
mean_time <- 1/mean_
text <- paste(mean_time,"days")
p <- ggplot(data=df,
aes(x = x,
y = gamma_s)) +
geom_line(colour = "#000000") +
geom_area(alpha = 0.6, fill = "lightgray") +
geom_vline(xintercept = mean_, linetype="dotted",
color = "red", size=0.5) +
annotate(geom="text", x=1.30 * mean_, y=0, label = text,
color="black")
return(p)
}
pi_mu_h <- mu_(lambda = 25, ylabel= "mu_h")
pi_mu_i_s <- mu_(lambda = 35, ylabel= "mu_i_s")
pi_gamma_s <- gamma_(alpha = 10, beta = 100, ylabel= "gamma_S")
pi_gamma_a <- gamma_(alpha = 10, beta = 50, ylabel= "gamma_A")
pi_gamma_h <- gamma_(alpha = 10, beta = 250, ylabel= "gamma_H")
fig6 <- ggplotly(pi_mu_h)## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
fig7 <- ggplotly(pi_mu_i_s)
fig8 <- ggplotly(pi_gamma_s)
fig9 <- ggplotly(pi_gamma_a)
fig10 <- ggplotly(pi_gamma_h)
fig6fig7fig8fig9fig10| theta | mean | se_mean | sd | 2.5% | 25% | 50% | 75% | 97.5% | n_eff | Rhat |
|---|---|---|---|---|---|---|---|---|---|---|
| beta_s | 8.690483e-01 | 0.0051487 | 0.0978857 | 6.821660e-01 | 8.046863e-01 | 8.683529e-01 | 9.370239e-01 | 1.061825e+00 | 361.44444 | 1.0014380 |
| beta_a | 7.738431e-01 | 0.0046027 | 0.1046186 | 5.690834e-01 | 7.016002e-01 | 7.722679e-01 | 8.464557e-01 | 9.786888e-01 | 516.63452 | 0.9995944 |
| kappa | 2.143749e-01 | 0.0041385 | 0.0453536 | 1.506806e-01 | 1.847692e-01 | 2.042473e-01 | 2.344915e-01 | 3.419625e-01 | 120.09807 | 1.0212692 |
| p | 4.387130e-01 | 0.0024290 | 0.0526256 | 3.089264e-01 | 4.087240e-01 | 4.528404e-01 | 4.803055e-01 | 4.989015e-01 | 469.38544 | 1.0011927 |
| delta_H | 4.603172e-01 | 0.0123912 | 0.1174493 | 1.971069e-01 | 3.909769e-01 | 4.616815e-01 | 5.355033e-01 | 6.810346e-01 | 89.84121 | 1.0140202 |
| mu_{I_S} | 2.602712e-01 | 0.2386425 | 1.3275922 | 4.162000e-03 | 2.803170e-02 | 5.691990e-02 | 1.042763e-01 | 3.605100e-01 | 30.94808 | 1.0442815 |
| mu_H | 2.327870e-01 | 0.2376622 | 1.3076330 | 2.533600e-03 | 1.305290e-02 | 3.280390e-02 | 6.643520e-02 | 3.525362e-01 | 30.27273 | 1.0449172 |
| gamma_S | 1.227483e-01 | 0.0034395 | 0.0412130 | 5.379330e-02 | 9.309530e-02 | 1.197828e-01 | 1.465723e-01 | 2.144921e-01 | 143.57880 | 1.0148258 |
| gamma_A | 5.106448e-01 | 0.0079512 | 0.0939779 | 3.163330e-01 | 4.584978e-01 | 5.113880e-01 | 5.662136e-01 | 6.887404e-01 | 139.69809 | 1.0283537 |
| gamma_H | 5.079869e-01 | 0.0089289 | 0.1639180 | 2.524056e-01 | 3.875971e-01 | 4.880279e-01 | 6.033939e-01 | 8.925461e-01 | 337.02411 | 1.0089564 |
| R_zero | 1.916769e+00 | 0.0988256 | 0.6028527 | 1.434153e+00 | 1.656316e+00 | 1.797068e+00 | 1.978445e+00 | 3.277463e+00 | 37.21206 | 1.0369071 |
| S_0 | 8.916145e+06 | 66.4756579 | 606.5895810 | 8.914838e+06 | 8.915746e+06 | 8.916239e+06 | 8.916569e+06 | 8.917209e+06 | 83.26541 | 1.0080175 |
| E_0 | 1.686677e+03 | 15.9085126 | 326.3754450 | 8.965301e+02 | 1.486263e+03 | 1.766709e+03 | 1.949267e+03 | 2.084575e+03 | 420.89697 | 1.0034477 |
| I_{S_0} | 2.200000e+01 | NA | 0.0000000 | 2.200000e+01 | 2.200000e+01 | 2.200000e+01 | 2.200000e+01 | 2.200000e+01 | NA | NA |
| I_{A_0} | 7.993784e+02 | 72.3470929 | 554.4589875 | 1.056505e+02 | 3.532306e+02 | 6.434443e+02 | 1.145775e+03 | 2.035521e+03 | 58.73496 | 1.0048948 |
| H_0 | 0.000000e+00 | NA | 0.0000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | NA | NA |
| R_0 | 0.000000e+00 | NA | 0.0000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | NA | NA |
| D_0 | 0.000000e+00 | NA | 0.0000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | NA | NA |
| ICS_0 | 7.400000e+01 | NA | 0.0000000 | 7.400000e+01 | 7.400000e+01 | 7.400000e+01 | 7.400000e+01 | 7.400000e+01 | NA | NA |
\[\begin{align} S'(t) &= \mu \bar{N} -\frac{\beta_S I_S+\beta_AI_A}{\bar{N}}S -(\mu+\lambda_V)S +\delta_V V \\ E'(t) &= \frac{\beta_S I_S + \beta_A I_A}{\bar{N}} S + \epsilon \frac{\beta_S I_S + \beta_A I_A}{\bar{N}}V - (\mu+\delta_E) E \nonumber \\ I'_S(t) & = p \delta_E E -(\mu + \mu_S + \alpha_S + \lambda_T) I_S +(1 - q) \alpha_T T \\ I'_A(t) & = (1-p) \delta_E E - (\mu + \mu_A + \alpha_A) I_A \\ R'(t) & = \alpha_S I_S + \alpha_A I_A + q \alpha_T T - \mu R \\ D'(t)& = \mu_S I_S + \mu_A I_A \\ V'(t)& = \lambda_V S - \epsilon \frac{\beta_S I_S + \beta_A I_A}{\bar{N}} V - (\mu + \delta_V) V \\ T'(t) &= \lambda_T I_S -(\mu + \alpha_T) T \\ \bar{N}(t) &= S(t) + E(t) + I_S(t) + I_A(t) + R(t) + V(t) + T(t) \end{align}\]
| state | Stan | Initial Condition | Stan | Postulated |
|---|---|---|---|---|
| \(S\) | y[1] | \(S_0\) | y_init[1] | \(N - (I_{S_0} + I_{A_0} + E_0)\) |
| \(E\) | y[2] | \(E_0\) | y_init[2] | Unif(0, 10) |
| \(I_s\) | y[3] | \(I_{s_0}\) | y_init[3] | Unif(0, 3) |
| \(I_a\) | y[4] | \(I_{a_0}\) | y_init[4] | Unif(0, 10) |
| \(R\) | y[5] | \(R_0\) | y_init[5] | Estimated |
| \(T\) | y[5] | \(R_0\) | y_init[5] | 0 |
| \(V\) | y[5] | \(R_0\) | y_init[5] | 0 |
| \(I_s^{[cumulative]}\) | y[6] | \(I_{s_0}^{[cumulative]}\) | y_init[6] | \(I_{s_0}\) |
| stan notation | symbol | Prior | Fixed |
|---|---|---|---|
| theta[1] | \(\beta_A\) | Normal(1, 0.3) | |
| theta[2] | \(\beta_S\) | Normal(1, 0.3) | |
| theta[3] | \(p\) | Unif(.3, 0.8) | |
| theta[3] | \(q\) | Unif(0.4, 0.7) | |
| theta[4] | \(\epsilon\) | 0.0 | |
| theta[5] | \(\alpha_S\) | Gamma(10, 100) | |
| theta[6] | \(\alpha_A\) | Gamma(10, 50) | |
| theta[7] | \(\alpha_T\) | 0.0 | |
| theta[8] | \(\lambda_V\) | 0.0 | |
| theta[9] | \(\lambda_T\) | 0.0 | |
| theta[10] | \(\mu_S\) | ||
| ——— | \(\mu_A\) | 0.0 | |
| ——— | \(\mu\) | ——– | 3.913894e-05 |
| theta[11] | \(\delta_V\) | 0.0 | |
| theta[12] | \(\delta_E\) | Gamma(10, 50) | 0.1960784 |
| ——— | \(E_0\) | Unif(0, 10) | |
| ——— | \(I_A\) | Unif(0, 10) | |
| ——— | \(I_S\) | Unif(0, 10) | |
| ——— | \(I_A\) | Unif(0, 10) |
| variable | Stan | Initial Condition | Stan | |
|---|---|---|---|---|
| \(S\) | y[1] | \(S_0\) | y_init[1] | \(N - (I{s_0} + I_{a_0} + E_0)\) |
| \(E\) | y[2] | \(E_0\) | y_init[2] | Unif(0, 10) |
| \(I_s\) | y[3] | \(I_{s_0}\) | y_init[3] | Unif(0, 3) |
| \(I_a\) | y[4] | \(I_{a_0}\) | y_init[4] | Unif(0, 10) |
| \(R\) | y[5] | \(R_0\) | y_init[5] | 0 |
| \(D\) | y[6] | \(D_0\) | y_init[6] | 0 |
| \(V\) | y[7] | \(V_0\) | y_init[7] | 0 |
| \(T\) | y[8] | \(V_0\) | y_init[8] | 0 |
| \(I_s^{[cumulative]}\) | y[9] | \(I_{s_0}^{[cumulative]}\) | y_init[9] | \(I_{s_0}\) |
## functions {
## real[] seirvt(real t, // time
## real[] y,
## real[] theta,
## real[] x_r,
## int[] x_i) {
## //
## real dy_dt[8];
## real s = y[1];
## real e = y[2];
## real i_s = y[3];
## real i_a = y[4];
## real h = y[5];
## real r = y[6];
## real d = y[7];
## real n_star;
## //
## real beta_s = theta[1];
## real beta_a = theta[2];
## real kappa = theta[3];
## real p = theta[4];
## real delta_h = theta[5];
## real mu_i_s = theta[6];
## real mu_h = theta[7];
## real gamma_s = theta[8];
## real gamma_a = theta[9];
## real gamma_h = theta[10];
## //
## real force_infection;
## real mu = 3.913894e-05;
## n_star = s + e + i_s + i_a + h + r;
##
## force_infection = (beta_s * i_s + beta_a * i_a) / n_star;
## dy_dt[1] = mu * n_star - force_infection * s - mu * s;
## dy_dt[2] = force_infection * s - (mu + kappa) * e;
## dy_dt[3] = p * kappa * e - (gamma_s + mu_i_s + delta_h) * i_s;
## dy_dt[4] = (1.0 - p) * kappa * e - (gamma_a + mu) * i_a;
## dy_dt[5] = delta_h * i_s - (gamma_h + mu_h + mu) * h;
## dy_dt[6] = gamma_s * i_s + gamma_a * i_a + gamma_h * h - mu * r;
## dy_dt[7] = mu_i_s * i_s + mu_h * h;
## dy_dt[8] = p * kappa * i_s;
## return dy_dt;
## }
## }
## data {
## int<lower = 1> n_obs; // number of days observed
## int<lower = 1> n_theta; // number of model parameters
## int<lower = 1> n_difeq; // number of differential equations
## int<lower = 1> n_pop; // population
## int y[n_obs]; // data, total number of infected
## real t0; // initial time point (zero)
## real ts[n_obs]; // time points observed
## }
##
## transformed data {
## real x_r[0];
## int x_i[0];
## }
## parameters {
## real <lower = 5e-6> theta[n_theta];
## // real <lower = 0, upper = 1> E0;
## // real <lower = 0, upper=1> IA0;
## real <lower = 0, upper = n_pop> E0;
## real <lower = 0, upper = n_pop> IA0;
## }
## transformed parameters{
## real y_hat[n_obs, n_difeq]; // solution from the ODE solver
## real y_init[n_difeq];
## // initial conditions
## // real IS0 = 8.297217e-06;
## real IS0 = 22.0;
## real H0 = 0.0;
## real R0 = 0.0;
## real D0 = 0.0;
## real CIS0 = 74.0;
## y_init[1] = n_pop - (IS0 + IA0 + E0);
## y_init[2] = E0;
## y_init[3] = IS0;
## y_init[4] = IA0;
## y_init[5] = H0;
## y_init[6] = R0;
## y_init[7] = D0;
## y_init[8] = CIS0;
## y_hat = integrate_ode_rk45(seirvt, y_init, t0, ts, theta, x_r, x_i);
## }
## model {
## real lambda[n_obs]; // poisson parameter
## // priors
## theta[1] ~ normal(1.0, 0.1); // beta_s
## theta[2] ~ normal(1.0, 0.1); // beta_a
## theta[3] ~ gamma(10, 40); // kappa
## theta[4] ~ uniform(0.1, 0.5); // p
## theta[5] ~ gamma(10, 40); // delta_h
## theta[6] ~ exponential(33.33333); // mu_i_s
## theta[7] ~ exponential(25); // mu_h
## theta[8] ~ gamma(10, 100); // gamma_s
## theta[9] ~ gamma(10, 50); // gamma_a
## theta[10] ~ gamma(10, 250); // gamma_h
## //
## // E0 ~ uniform(1.121246e-07, 3.363737e-05);
## // IA0 ~ uniform(7.848719e-06, 3.363737e-05);
## E0 ~ uniform(74, 2100);
## IA0 ~ uniform(74, 2100);
## // likelihood
## for (i in 1:n_obs){
## //lambda[i] = y_hat[i, 8] * n_pop;
## lambda[i] = y_hat[i, 8];
## }
## y ~ poisson(lambda);
## }
## generated quantities {
## real R_0; // Basic reproduction number
## real mu = 3.913894e-05;
## real beta_s = theta[1];
## real beta_a = theta[2];
## real p = theta[4];
## real gamma_s = theta[4];
## real gamma_a = theta[5];
## real kappa = theta[6];
## R_0 = p * beta_s * kappa / ((gamma_s + mu) * (kappa + mu))
## + (1.0 - p) * beta_a * kappa / ((gamma_a + mu) * (kappa + mu));
## }
\[\begin{align} Y_t & \sim \text{Poisson}(\lambda_t), \qquad \\ \lambda_t & = \int_0^t p \kappa E \\ & p \sim \text{Uniform(0.3 0.8)} \\ & \kappa \sim \text{Gamma(10, 52)} \end{align}\]